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1.  INTRODUCTION 


To  determine  the  accuracy  of  artillery  fire  one  measures  the  coordinates  of  the 
shell’s  burst  point  in  repeated  firings  and  calculates  an  average  burst  point  and  its 
scatter  from  these  measurements.  The  task  amounts  to  the  computation  of  an  average 
vector  in  The  accuracy  of  each  observed  vector  is  known  from  an  analysis  of  the 
actual  measurements  and  depends  mainly  on  the  geometry  of  the  setup  and  properties  of 
the  measuring  instruments  (theodolites  in  general).  We  characterize  this  accuracy  by  an 
estimated  variance-covariance  matrix  of  the  vector  components.  If  the  cannon  would 
fire  every  time  exactly  alike  (i.e.,  if  the  event  scatter  would  be  zero)  then  the  burst- point 
coordinates  could  be  obtained  from  these  data  by  a  weighted  averaging  where  the 
weights  are  the  inverses  of  the  variance-covariance  matrices.  However,  in  real  life  the 
event  scatter  can  be  of  the  same  order  or  even  larger  than  the  measurement  scatters. 
Also,  in  general  the  principal  directions  of  the  event  distribution  are  different  from  the 
principal  directions  of  the  measurement-error  distributions.  Therefore,  a  weighted 
averaging  in  can  produce  unacceptable  results.  On  the  other  hand  an  unweighted 
averaging  would  not  take  into  account  the  estimates  of  measurement  errors  that  can  be 
quite  different  for  different  observations.  In  this  paper  we  present  a  new  algorithm  for 
the  computation  of  an  average  vector  that  does  not  have  the  disadvantages  of 
unweighted  or  observation-weighted  averages.  The  algorithm  provides  in  addition  to 
the  average  vector,  also  an  estimate  of  the  event  variance  that  is  consistent  with  the 
observations  and  their  estimated  variances. 

In  Section  2  we  define  a  problem  of  vector  averaging  in  SV*  that  corresponds  to  the 
outlined  artillery  problem  and  propose  a  solution.  Section  3  contains  some  examples 
and  Section  4  is  a  summary. 

2.  ESTIMATION  OF  AN  AVERAGE  VECTOR 

Let  the  observed  vectors  be  z,-  €  f?",  i  =  1,...,5  and  let  the  estimated  variance- 
covariance  n  X  n  matrices  of  the  observations  be  Q^,  i  =  Let  the  unknown 

average  vector  be  a  6  JF?"  and  the  variance-covariance  matrix  of  the  event  be  P.  The 
model  equation  of  the  event  is 

/(z,o)  =  z  -  0  =0  .  (1) 

We  define  the  least-squares  value  of  a  as  the  solution  of  the  following  constrained 
optimization  problem. 
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Minimize 


(2) 


yV^ticTQr\  +  bTp-%) 

1—1 

subject  to 

/(*,•  +  c,-.  a  +  6.)  =  I, •  +  c,  -  a  -  fc,-  =0  ,  i  =  , 

where  c,-  is  the  correction  of  the  i-th  observation  and  6,-  is  the  deviation  of  the  i-th  event 
from  the  average  a. 

To  solve  the  minimization  problem  we  introduce  a  modified  objective  function  W 
using  Lagrange  multiplier  vectors  A;,-: 


S  {c^Qr^Ci  +  bJP  )-Yi  kj  {Xi  +  Ci 

1-1  »-i 


a-bi)  .  (4) 


We  obtain  a  system  of  normal  equations  by  setting  equal  to  zero  the  partial  derivatives 
of  W  with  respect  to  c^,  6,-,  a  and  k^.  The  result  is 


Qr'ci  - 1. 

p-H,  +  ki 


Xi  +  Ci  —  a  -  bi 


i  = 


(5) 


Eliminating  the  we  obtain  the  following  simpler  equation  system 

a  =  [  E  (g.-  +  P)-^  E  {Qi  +  Pp  X.-  . 

6,.  =  P{Qi+Pr^{xi-a)  ,  i  =  l . s, 

Ci  =-  Qii  Qi-^  P)~^  (xi-  a)  ,  i  =  . 


(6) 


We  also  obtain 


w,  =  t  iJP-'K  =  t  (*.-  -  «)^(  <2,  +  /■)-■/’  (  Qi  +  P )-'  (I.-  -  0)  ,  (7) 

i— 1  1—1 

K  =  t  cTq-\  =  t  (%  -  «)^(  Qi  +  P  )■'  Qi  (Qi  +  P  )■'  (*.•  -  “)  .  (*) 

•-1  »-l 

1—1 


and  the  variance  of  weight  one 


-2- 


(10) 


‘'o  —  /  1  \  • 

n  («  -  1) 

Let  the  total  variance-covariance  matrix  including  both  the  measurement  scatter 
and  the  event  scatter  of  the  observed  z,-  be  Then  the  variance-covariance  matrix  of 


«.=[£(%+/>)-']  ‘  4  [  (<3, + pr'R.i  (Q. + f)-'  I  lijQi+p  r'l"’-  (u) 

i—i  1—1  1—1 


If  we  estimate  as  usual 

then  it  follows  from  eq.  (11) 


The  formulas  (6)  through  (13)  provide  the  general  least-squares  solution  of  the  averaging 
problem  if  the  and  P  are  known.  In  practice  such  a  situation  is  an  exception,  because 
ip  general  P  is  not  known.  Therefore,  commonly  used  are  two  special  cases  of  the 
solution  that  are  based  on  the  assumption  that  either  the  Q,-  or  P  can  be  neglected.  We 
now  outline  these  special  cases. 

In  the  first  special  case  we  assume  that  P  =  0,  that  is,  we  assume  that  either  the 
event  scatter  is  negligible  or  that  the  estimated  Q,-  already  contain  the  matrix  P.  With 
this  assumption  we  obtain  from  eqs.  (6)  and  (13)  the  usual  observation- weigh  ted  least- 
squares  averaging  formulas: 

a  =  [E  Q.-'l"'  E  QrS  .  (14) 

hj  =0  ,  I  =  1,...,5  , 

c,  =  a  -  ij  ,  i  =  !,...,«  , 

«.  =  »„[  <?r' )''  ■  (>6) 

Usually  the  Q,-  are  positive  definite  matrices  but  in  some  applications  they  may  be  only 
semi-definite.  Also,  the  sum  of  their  inverses  in  eqs.  (14)  and  (16)  is  not  necessarily 
positive  definite.  The  formulas  are,  however,  generally  valid  if  Moore-Penrose 
generalized  inverses  are  used  in  both  formul 

In  the  second  case  we  assume  that  the  measurement  errors  are  negligible,  that  is, 
Qi—0  for  all  i  =  !,...,«  (or  that  all  Q,-  are  equal  and  included  in  P).  In  that  case 
eqs.  (6)  and  (13)  provide  the  formulas 
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(17) 


1  ^ 

S  1-1 

bi  =  Xi-a  , 

c,  =0  , 


1=1,.. .,5  , 
j  =  , 


«.  =  «o  7  ^  •  (19) 

Equations  (17)  through  (19)  are  the  formulas  for  simple  unweighted  averaging.  To 
complete  the  calculation  we  also  need  an  estimate  for  P.  The  usual  estimate  is 

P=aP  =  aEbibl  ,  (20) 

where  the  factor  a  is  determined  such  that  equals  unity.  To  that  end  we  compute 

^=11  (i,  -  af  P-^  (Xi  -  a)  ,  (21) 

»— 1 

^  ^  U  (22) 

n  (s  -  1) 

and  obtain  for  the  variance-covariance  matrix  of  the  average 

R,  =  v,^P - (23) 

3  n  3  (3  —  1) 

As  in  the  first  case,  we  use  the  Moore-Penrose  generalized  inverse  in  eq.  (21)  if  the 
matrix  P  is  not  positive  definite. 

In  real  life,  estimates  of  the  Q,-  usually  are  available  but  P  is  not  known  so  that 
the  general  formulas  cannot  be  used.  If  also  neither  of  the  two  special  cases  apply  one 
needs  a  method  that  provides  an  estimate  of  P  concurrently  with  the  average  vector  o. 
We  propose  for  this  purpose  the  following  iterative  algorithm.  It  produces  an  estimate 
of  P  and  solves  the  general  problem,  eqs.  (2)  and  (3),  using  this  estimate.  Because  the 
solution  takes  into  account  the  distribution  of  the  observed  vectors  z,-  we  call  the 
resulting  a  the  distribution-weighted  average. 

We  initialize  the  iteration  with  an  unweighted  averaging 


b^i^Xf-a,  1  =  1,. ..,3, 


and  obtain  an  initial  estimate  of  P 


-4- 


(26) 


s 


'0,1  ''0,1 


Uo  =  g  bl,  fc,.,  , 

F  P 

nis-l)  ■ 


(27) 

(28) 


Next,  we  update  the  initial  estimates  of  a  and  6,-  and  obtain  an  initial  estimate  for  the 
scaling  factor  a. 


<■1 = 1 1:  (ft + p,r'  1  ‘  g  (ft + ^i)-'  *. .  (29) 

‘m  -  ^>1  (ft +  f’i)'’ (»/-<■.),  .  . .  (30) 

P I  “  g  4f  j  ,  (31) 

.  (92) 

n  («  -  1) 

In  subsequent  iteration  steps  we  do  update  P  but  keep  the  value  of  or  unchanged.  The 
iteration  formulas  for  k  3=  1,2,...  are 

^*+1  “ 

0,«  =  [  t  (ft  +  Pm)-'  r‘  E  (ft  +  PmT'  .  (35) 

h+I.$  ^  Pt+l  Pt*l)  *  (*i  “  ®t+l)  .  ’  =  1 . ®  >  (36) 

Pm  =  t  iM.i  ‘.♦I,.-  ■  (37) 

•—1 


The  variance-covariance  estimate  F,  of  the  average  is  computed  with  eqs.  (9),  (10)  and 
(13).  Iteration  end  conditions  can  be  expressed,  for  instance,  in  terms  of  changes  of  the 
elements  of  a  and  F,.  Experience  shows  that  the  average  vector  a  becomes  stationary 
after  a  few  iteration  steps  whereas  the  elements  of  R,  need  more  steps  to  meet  such 
convergence  criteria.  Convergence  enhancement  techniques  were  found  unnecessary  in 
numerical  experiments  with  this  algorithm. 


The  result  of  the  iterations  depends  on  the  scaling  factor  a  that  was  initially 
estimated  by  eq.  (33).  We  want  to  determine  its  value  such  that  the  variance  of  weight 
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one  Vg,  defined  by  eqs.  (9)  and  (10),  equals  unity.  We  achieve  this  by  embedding  the 
iteration  eqs.  (34)  through  (37)  in  a  regula  falsi  algorithm  for  the  solution  of  the 
equation  Vo(®)  =  1-  1“  general,  a  solution  with  positive  a  exists  if  ^^(0)  >  1,  because 
decreases  with  increasing  a.  If  vjO)  <  1  then  the  estimated  observational  errors  (the 
matrices  Q,)  are  so  large  that  a  solution  with  or  0,  i.e.,  with  neglected  P  suffices  to 
explain  the  data. 

The  final  result  of  the  calculations  is  a  solution  of  the  general  minimization 
problem,  eqs.  (2)  and  (3),  whereby  the  event  variance  matrix  satisfies  eq.  (20)  and 
v„  =  l. 

3.  EXAMPLES 

We  present  two  examples.  The  first  example  is  chosen  to  illustrate  the  main 
characteristics  of  the  three  types  of  averaging.  In  the  second  example  we  use  actual 
data. 


In  the  first  example  we  compute  the  average  of  three  points  on  a  straight  line  in  a 
plane.  The  coordinates  of  the  points  are  (0.5,  2.0),  (1.5,  2.1)  and  (8.5,  2.8).  We  assume 
that  the  observational  errors  are  equal  for  all  three  points  and  given  by  the  following 
estimate  of  their  variance-covariance  matrix 


2.0  ) 

2.0  J  • 


The  matrix  Q  is  not  positive  definite  which  means  that  the  observational  errors  are 
distributed  in  a  subspace  of  IR^,  that  is,  along  a  straight  line.  In  other  words,  the 
observational-error  ellipses  are  degenerated  into  error  bars.  Figure  1  shows  the  data 
and  the  result  of  a  weighted  average.  The  coordinates  of  the  average  are  (3.5,  2.3)  and 
the  variance-covariance  matrix  of  the  average  is 

P  _  {  0.95792  0.95792  ) 

“  I  0.95792  0.95792  J  ' 

The  corresponding  standard-deviation  ellipse  is  again  degenerated  and  shown  in  Figure  1 
as  a  dashed  line.  The  location  of  the  average  point  is  reasonable  but  its  estimated 
variance  is  not  because  the  structure  of  the  variance-covariance  matrix  that  is  computed 
with  eq.  (16)  is  independent  of  the  observations  and  does  not  reflect  the  event  scatter. 


Next  we  use  the  same  data  and  compute  their  unweighted  average  by  eqs.  (17) 
through  (23).  The  average  vector  is  the  same  as  in  the  previous  calculation  but  its 
variance-covariance  matrix  is 


p  _  (  4.75  0.475  ) 

•  “  I  0.475  0.0475  J  ' 

The  result  is  shown  in  Figure  2.  The  image  of  the  one-standard-error  ellipse  of  the 
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Figure  1.  Observation- weigh  ted  average. 


average  is  an  error  bar  in  the  direction  of  the  scatter  of  the  observations,  because  in  this 
case  is  independent  of  the  observational-error  variances. 


Figure  2.  Unweighted  average. 

Finally,  we  compute  the  distribution-weighted  average  by  iteration,  eqs.  (24) 
through  (37).  The  average  vector  again  is  the  same  as  before.  Its  variance-covariance 
matrix  is 


-7- 


p  _  f  4.12682  1.13567  ] 

•  ~  [  1.13567  0.73027  J  ' 

Figure  3  shows  the  corresponding  one-standard-error  ellipse.  The  figure  also  contains 
the  correction  v?  -tors  6,-  plotted  as  rays  from  the  average  point.  The  end  points  of  the 
6,-  are  indicated  by  dots.  In  this  example,  all  vectors  5,-  are  parallel  so  that  their  end 
points  are  located  on  a  straight  line  and  the  matrix  P,  eq.  (34),  is  only  positive  semi- 
definite.  The  image  of  the  ellipse  representing  P  is  a  segment  of  a  straight  line  in  the 
direction  of  the  6,-.  The  differences  between  the  dots  and  the  corresponding  observations 
(inverted  triangles)  are  the  corrections  c^.  We  observe  that  all  corrections  c,-  and  6,-  are 
in  the  direction  of  the  corresponding  error  bars,  as  they  should  be.  In  this  example  the 
iteration  with  eqs.  (34)  through  (37)  became  stationary  after  two  steps.  The  initial 
scaling  factor  and  the  variance  of  weight  one  were,  respectively,  a  =  0.250  and 
Uq  =  1.008.  After  three  regula  falsi  steps,  we  had  the  values  a  =0.252006  and 

Vo  =  1.00006. 


East,  m 

Figure  3.  Distribution-weighted  average. 

In  our  second  example  we  use  actual  observations  of  artillery  burst-point 
coordinates.  The  observations  are  three-dimensional  vectors  consisting  of  range, 
deflection  and  height  of  the  burst.  The  coordinates  of  each  vector  were  obtained  from 
simultaneous  measurements  of  directional  angles  of  the  burst  points  from  four 
observation  towers.  An  analysis  of  these  measurements  provided  the  components  of 
each  burst  location  vector  and  estimates  of  the  accuracies  of  the  burst  points  in  form  of 
the  variance-covariance  matrices  Qy  The  estimated  accuracies  of  the  observations  vary 
widely  and  are  in  this  example  smaller  than  the  scatter  between  the  burst  points,  but 
are  not  negligible.  The  observation  set  in  our  example  consists  of  eight  observed  burst 
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points  from  the  same  howitzer.  Figures  4,  5  and  6  show  the  observed  points  as  inverted 
triangles  and  their  distribution-weighted  (iteratively  determined)  average  as  a  diamond. 
The  figures  also  contain  the  projections  of  the  corresponding  one-standard-error 
ellipsoids.  The  standard-error  ellipsoid  of  the  average,  defined  by  R^,  is  plotted  with  a 
dashed  line.  The  standard-deviation  ellipsoid  of  a  single  shot  is  represented  by  the 
matrix  P  and  plotted  with  a  dotted  line.  We  note  that  contrary  to  the  appearance  in 
the  plots  P  is  not  proportional  to  R,:  the  relation  between  P  and  R,  is  given  by 
eq.  (13).  The  correction  vectors  6,-  are  plotted  as  rays  from  the  average  point,  as  in 
Figure  3.  We  observe  that  these  corrections  are  in  general  not  in  the  direction  towards 
the  observations,  but  towards  other  points  such  that  the  corrections  6,-  and  c,-  are  in 
directions  of  large  error  estimates  thus  minimizing  W,  eq.  (2).  The  initial  value  of  the 
scaling  factor  was  0=0.143  for  a  variance  of  weight  one  of  =  1.137.  After  four 
regula  falsi  steps,  the  results  were  ot  =0.176998  and  v,,  =  1.00004.  The  iteration  for  a 
and  6,-,  eqs.  (34)  through  (37),  required  eight  steps  at  the  beginning  and  three  steps  at 
the  end  of  the  calculations. 


Firing  range,  m 

Figure  4.  Burst-point  range  and  deflection. 

To  illustrate  the  advantage  of  the  distribution-weighted  average  we  show  in 
Figures  7,  8  and  9  the  usual  observation-weighted  average  [eqs.  (14)  through  (16)  ]  of  the 
same  observations.  We  notice  that  in  this  example  the  result  is  completely  unrealistic 
because  the  observation-weighted  average  burst  point  is  shifted  far  outside  the  cloud  of 
observations.  From  an  inspection  of  the  figures,  we  conclude  that  this  shift  is  caused  by 
the  high  sensitivity  of  the  location  of  the  average  to  the  estimated  principal  directions  of 
observational  errors.  The  variance  of  weight  one  was  in  this  case  =  5996  indicating 
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Figure  5.  Burst- point  range  and  height. 


Figure  6.  Burst-point  deflection  and  height, 
that  measurement  errors  alone  are  not  sufficient  to  explain  the  data  scatter. 

4.  SUMMARY 

We  have  considered  least-squares  computations  of  vector  averages.  We  assume 
that  the  observations  in  an  n-dimensional  space  contain  inaccuracies  from  two  sources; 
observational  errors  and  variations  of  the  observable  itself,  that  is,  event  scatter. 
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Figure  7.  Observation-weighted  range  and  deflection. 
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Figure  8.  Observation-weighted  range  and  height. 

Usually  one  of  these  error  sources  is  neglected.  If  a  simple  unweighted  average  is 
computed  then  one  assumes  implicitly  that  the  observational  errors  are  negligible.  If  a 
weighted  average  is  computed  then  the  implicit  assumption  is  that  the  event  scatter  can 
be  neglected.  Most  often  the  event  variation  is  not  known  and  one  has  no  grounds  for 
using  one  of  these  special  algorithms.  If  event  scatter  is  known  to  exist  we  propose  to 
use  the  distribution- weighted  average  that  is  computed  by  an  iterative  algorithm.  The 


_ _ 

\ 

f 

(' 

N 

< 

_ 

-  11  - 
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Figure  9.  Observation- weighted  deBection  and  height. 

algorithm  provides  in  addition  to  the  average  of  the  observed  vectors  with  its  variance, 
also  an  estimate  of  the  variance-covariance  matrix  of  the  event. 
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